clear all;
clc;

% ---------------------- Import Data --------------------------------------

cd('patient1_8_12_13');
dc_xy=load('p1_dc_xyzloc.txt');
dr_xy=load('p1_dr_xyzloc.txt');


dc_xyz=load('p1_dcolmxyz_200axons_312nodes_9um.txt');
dr_xyz=load('p1_drootxyz_200axons_167nodes_9um.txt');


% cd('other_geom');
% atpphi_dc = load('sub_atp_dc_phi_p1_d0p2_t0.txt');
% atpphi_dr = load('sub_atp_dr_phi_p1_d0p2_t0.txt');
% cd ..;

cd('longtripole_6ies');
atpphi_dc = load('sub_LT6_dc_phi_p1_d0p2_t0.txt');
atpphi_dr = load('sub_LT6_dr_phi_p1_d0p2_t0.txt');
cd ..;

cd ..;

% ---------------------- Process Data -------------------------------------

[noaxons_dc,nopts_dc]=size(atpphi_dc);

dc_xc=zeros(nopts_dc,noaxons_dc);
dc_yc=zeros(nopts_dc,noaxons_dc);
dc_zc=zeros(nopts_dc,noaxons_dc);
dc_xc(:)=dc_xyz(1,:);
dc_yc(:)=dc_xyz(2,:);
dc_zc(:)=dc_xyz(3,:);

[noaxons_dr,nopts_dr]=size(atpphi_dr);

dr_xc=zeros(nopts_dr,noaxons_dr);
dr_yc=zeros(nopts_dr,noaxons_dr);
dr_zc=zeros(nopts_dr,noaxons_dr);
dr_xc(:)=dr_xyz(1,:);
dr_yc(:)=dr_xyz(2,:);
dr_zc(:)=dr_xyz(3,:);

% analyze only nodes of Ranvier
phi_dcnodes=atpphi_dc(:,1:8:end);
sd_dcnodes=diff(phi_dcnodes,2,2);
dcnodes_xc=dc_xc(1:8:end,:)';
dcnodes_yc=dc_yc(1:8:end,:)';
dcnodes_zc=dc_zc(1:8:end,:)';

dxtmp=dcnodes_xc(:,2:end)-dcnodes_xc(:,1:end-1);
dytmp=dcnodes_yc(:,2:end)-dcnodes_yc(:,1:end-1);
dztmp=dcnodes_zc(:,2:end)-dcnodes_zc(:,1:end-1);
dcnodes_arc=[zeros(noaxons_dc,1),cumsum(sqrt(dxtmp.^2+dytmp.^2+dztmp.^2),2)];

phi_drnodes=atpphi_dr(:,1:8:end);
sd_drnodes=diff(phi_drnodes,2,2);
drnodes_xc=dr_xc(1:8:end,:)';
drnodes_yc=dr_yc(1:8:end,:)';
drnodes_zc=dr_zc(1:8:end,:)';

dxtmp=drnodes_xc(:,2:end)-drnodes_xc(:,1:end-1);
dytmp=drnodes_yc(:,2:end)-drnodes_yc(:,1:end-1);
dztmp=drnodes_zc(:,2:end)-drnodes_zc(:,1:end-1);
drnodes_arc=[zeros(noaxons_dr,1),cumsum(sqrt(dxtmp.^2+dytmp.^2+dztmp.^2),2)];

% max second difference
maxsd_dc=max(sd_dcnodes,[],2);
maxsd_dr=max(sd_drnodes,[],2);

[junk1,dcsd_si]=sort(maxsd_dc);
[junk2,drsd_si]=sort(maxsd_dr);


% ---------------------- Selecting Examples -------------------------------

subdc=130:180;
subdr=1:50;

aa=30;

ifar_dc=dcsd_si(aa);
ifar_dr=drsd_si(aa);
sdfar_dc=sd_dcnodes(ifar_dc,:);
sdfar_dr=sd_drnodes(ifar_dr,:);

inear_dc=dcsd_si(end-aa);
inear_dr=drsd_si(end-aa);
sdnear_dc=sd_dcnodes(inear_dc,:);
sdnear_dr=sd_drnodes(inear_dr,:);

a_dcnear=dcnodes_arc(inear_dc,2:end-1);
a_drnear=drnodes_arc(inear_dr,2:end-1);
a_dcfar=dcnodes_arc(ifar_dc,2:end-1);
a_drfar=drnodes_arc(ifar_dr,2:end-1);

% ----------------------------- Plotting ----------------------------------

figure;
plot(maxsd_dc(dcsd_si),'k-','LineWidth',4);hold on;
plot(maxsd_dr(drsd_si),'k--','LineWidth',4);hold off;
ylabel('Maximum \Delta^2\Phi (V/mm)','FontSize',36);
xlabel('Element Number','FontSize',36);
set(gca,'FontSize',36,'YScale','Log','XTick',25:25:200);
legend('DC fibers','DR fibers','Location','N');

figure;
subplot(2,2,1),plot(a_dcnear(subdc),sdnear_dc(subdc),'k-',...
    'LineWidth',4);
title('DC','FontSize',36);
ylabel('\Delta^2\Phi (V/mm)','FontSize',36);
set(gca,'FontSize',30);
xlim([a_dcnear(subdc(1)),a_dcnear(subdc(end))]);
% ylim([-0.2,0.2]);
subplot(2,2,2),plot(a_drnear(subdr),sdnear_dr(subdr),'k--',...
    'LineWidth',4);
title('DR','FontSize',36);
set(gca,'FontSize',30);
xlim([a_drnear(subdr(1)),a_drnear(subdr(end))]);
% ylim([-0.02,0.02]);

subplot(2,2,3),plot(a_dcfar(subdc),sdfar_dc(subdc),'k-',...
    'LineWidth',4);
xlabel('Arc Length (mm)','FontSize',36);
ylabel('\Delta^2\Phi (V/mm)','FontSize',36);
set(gca,'FontSize',30);
xlim([a_dcfar(subdc(1)),a_dcfar(subdc(end))]);
% ylim([-2.1e-3,7e-3]);
subplot(2,2,4),plot(a_drfar(subdr),sdfar_dr(subdr),'k--',...
    'LineWidth',4);
xlabel('Arc Length (mm)','FontSize',36);
set(gca,'FontSize',30);
xlim([a_drfar(subdr(1)),a_drfar(subdr(end))]);
% ylim([-2.1e-3,7e-3]);